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We study several variants of the stochastic four-state rock-paper-scissors game or, equivalently, 
cyclic three-species predator-prey models with conserved total particle density, by means of Monte 
Carlo simulations on one- and two-dimensional lattices. Specifically, we investigate the infiuence of 
spatial variability of the reaction rates and site occupancy restrictions on the transient oscillations 
of the species densities and on spatial correlation functions in the quasi-stationary coexistence state. 
For small systems, we also numerically determine the dependence of typical extinction times on the 
number of lattice sites. In stark contrast with two-species stochastic Lotka-Volterra systems, we 
find that for our three-species models with cyclic competition quenched disorder in the reaction 
rates has very little effect on the dynamics and the long-time properties of the coexistence state. 
Similarly, we observe that site restriction only has a minor infiuence on the system's dynamical 
properties. Our results therefore demonstrate that the features of the spatial rock-paper-scissors 
system are remarkably robust with respect to model variations, and stochastic fiuctuations as well 
as spatial correlations play a comparatively minor role. 
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I. INTRODUCTION 



Understanding the origin of and maintaining biodiver- 
sity is of obvious paramount importance in ecology and 
biology In this context, paradigmatic schematic 

models of predator-prey interaction that build on the 
classic Lotka-Volterra system 0] have been widely 
studied. Specifically, systems with cyclic dominance of 
competing populations have been suggested to provide 
a mechanism to promote species diversity; there are 
also natural connections to evolutionary game theory fi- 
lial. A minimal yet non-trivial model for cyclic com- 
petition is the three-species cyclic predator-prey system 
with standard Lotka-Volterra predation interactions, es- 
sentially equivalent to the familiar rock-paper-scissors 
(RPS) game fj-fllt. This RPS system has, for exam- 
ple, been used to model the cyclic competitions between 
three subspecies of certain Californian lizards [l3, fl5| . 
and the coevolution of three strains of E. coli bacteria 
in microbial experiments |l6l. Other examples include 
coral reef invertebrates [17| and overgrowths by marine 
sessile organisms [IBIill- In this simple RPS model, one 
lets 'rock' (species A) smash 'scissors' (species B), 'scis- 
sors' cut 'paper' (species C), and 'paper' wrap 'rock'. 
Already for a non-spatial RPS system, the presence of in- 
trinsic stochastic fluctuations (reaction noise) makes the 
system eventually evolve to one of the three extinction 
states where only one species survives [23-!^. For ex- 
ample, if the reaction rates in the system are not equal. 
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one intriguingly observes the 'weakest' species, with the 
smallest predation rate, to survive, whereas the other 
two species always die out jH, [1^ ■ When the model is 
extended to include spatial degrees of freedom, say by 
allowing particles to hop to nearest-neighbor sites on a 
lattice and interact upon encounter, spatial fluctuations 
and correlations further complicate the picture. For in- 
stance, species extinction still prevails in one-dimensional 
RPS models [23 - [27j . but the system settles in a coexis- 
tence state when the species are efficiently mixed through 
particle exchange (but see also Ref. [231 )■ In contrast, 
two-dimensional RPS systems are characterized by co- 
existence of the competing species, and the emergence 
of complex spatio-temporal structures such as spiral pat- 
terns 22, 24, 27, 29-38]. Recently, Reichenbach et al. ex- 
tensively studied the four-state RPS model without con- 
servation law fs^-H^ , and it is now well-established that 
cycUc reactions in conjunction with diffusive spreading 
generate spiral patterns (when the system is sufficiently 
mixed). In model variants that incorporate conserva- 
tion of the total population density, on the other hand, 
spiral patterns do not occur [12, Hj, [ss'] ; also, when the 
species mobility is drastically enhanced through fast par- 
ticle exchange processes, the spiral patterns are destroyed 
as well, and the system eventually reaches an extinction 
state [33i. i38j|. 

This work is motivated by the following question: 
Which are the crucial model ingredients to be included 
in order to attain a further degree of realism? To this 
end, we carried out Monte Carlo simulation studies on 
the influence of the carrying capacity and environmen- 
tal inhomogeneity on the properties of a class of spatial 
RPS models where the total population size is conserved 
(zero-sum games) |, III [lollS,^!! Ill Hill, IMIl- 
In our stochastic lattice models, the carrying capacity 
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(i.e., the maximum population size that can be sustained 
by the environment) is implemented through site occu- 
pation number restrictions. Environmental variability is 
modeled through assigning local reaction rates that are 
treated as quenched random variables drawn from a uni- 
form distribution. Our extensive numerical study shows 
that carrying capacity and quenched disorder have little 
influence on the oscillatory dynamics, spatial correlation 
functions, and extinction times in the RPS model sys- 
tem. This demonstrates a quite remarkable robustness 
of this class of models. From a modeling perspective, 
this establishes the essential equivalence of rather dis- 
tinct model variants. We emphasize that this outcome 
is nontrivial, as is, for example, revealed by a compari- 
son with the two-species Lotka-Volterra system [10, |4l| , 
where spatially varying reaction rates may cause more 
localized clusters of activity and thereby enhance the fit- 
ness of both predator and prey species [13] ■ 

Naturally, lattice models should be viewed as coarse- 
grained representations of a metapopulation system 
where each lattice site or cell can be interpreted as a 
'patch' (or 'island') populated by a 'deme' (or 'local com- 
munity') [13, 13]. For the sake of simplicity (i.e. to try 
to understand the simplest possible systems before ven- 
turing further), we here restrict our presentation to RPS 
model variants that obey a conservation law (even though 
that has no particular ecological motivation). As empty 
lattice sites are allowed, we shall refer to our model sys- 
tem as a class of four-state RPS models with conservation 
law (if all sites are at most occupied by a single individ- 
ual, each of them can be in one of four states). While 
the presence or absence of the conservation of the total 
number of particles is crucial for the emergence of spiral 
waves in RPS systems |38|, this is not the case for the 
properties studied here. In fact, it turns out that our 
conclusions on the effects (or lack thereof) of limited car- 
rying capacity and random environmental influences on 
the transient population oscillations, spatial correlation 
functions, and species extinction times are common to 
models both with and without conservation laws [45| . 

Our paper is structured as follows: In Sec. II, we define 
our model, the stochastic four-state spatial rock-paper- 
scissors (RPS) game or cyclic three-species predator-prey 
system with conservation of total population density, and 
briefly review the results obtained from the mean-field 
rate equation approximation. In Sec. Ill, we introduce 
our Monte Carlo simulation algorithm and discuss the de- 
tailed model variants we have explored. We then present 
results for the species' time-dependent densities, asso- 
ciated frequency power spectra, and spatial correlation 
functions to analyze the influence of quenched spatial 
disorder in the reaction rates and site occupation re- 
striction on the temporal evolution and quasi-stationary 
states of this system, both in two dimensions and for a 
one-dimensional lattice. We also compare our numerical 
findings with the mean-field predictions, and obtain the 
mean extinction time (for the first species to die out) as 
function of system size. Finally, we provide a summary 



of our results and concluding remarks. 

II. MODEL AND RATE EQUATIONS 

The rock-paper-scissors (RPS) model describes the 
cyclic competition of three interacting species that we 
label A, B, and C. We consider the following (zero-sum 
0) predator-prey type interactions: 

A + B A + A with rate ka , 

B + C ^ B + B with rate h , (1) 

C + A C + C with rate ■ 

Note that these irreversible reactions strictly conserve 
the total number of particles. We remark that naturally 
other variants of the RPS dynamics could also be consid- 
ered; notably the four-state May-Leonard model which 
does not conserve the total particle density [46'| has at- 
tracted considerable attention, see e.g. Refs. [9, 34, 47] . 
As will be demonstrated elsewhere, the conclusions pre- 
sented here on the effects of carrying capacity and spatial 
reaction rate variability remain essentially unchanged for 
this system [4^ . To generalize the above reaction model 
to a spatially extended lattice version, we allow empty 
sites (as a fourth possible state) and let the reactions 
happen only between nearest neighbors. In addition, we 
introduce nearest-neighbor particle hopping with rate D 
(if at most one particle is allowed per lattice site, this pro- 
cess takes place only if an adjacent empty site becomes 
selected at each time step). 

Within the mean-field approximation, wherein any cor- 
relations and spatial variations are neglected, the follow- 
ing set of three coupled rate equations for homogeneous 
population densities a(t), b{t), and c(i), with fixed total 
population density a{t) + b{t) + c{t) = p — const describes 
the system's temporal evolution, 

dt ait) = a{t) [ka h{t) - fc, c{t)] , 

dth{t) = b{t)[kbc{t)-kaa{t)] , (2) 

dtcit) = c{t)[kca{t)-hb{t)] . 

These coupled rate equations possess a reactive fixed 
point, where all three species coexist, {a*,b*,c*) — 
(kh, kr., kn.)p/(kn + kb + kc), which is marginally stable 
(see also Ref. [23). Indeed, introducing new variables 
Sa{t) = a{t) - a*, 5b{t) = b[t) - b* , 5c{t) = c{t) - c* , and 
utilizing the conservation law 6a + Sb + 6c — 0, we may 
express the first two rate equations in terms of i5a and 
6b. Linearizing about the reactive fixed point then gives 

with the linear stability matrix 
J- ^ P ( h kc kb {ka + kc) \ , . 

ka + kb + kc\~kc{ka+kb) -kbkc J ' 
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with eigenvalues A = ±i p \/ka kb kc/{ka + kb + kc) = 
zbia;, where / = lojl-K represents a characteristic os- 
cillation frequency, e.g., for total density p = 1 and 
/ca = 0.2, kb = 0.5, fee = 0.8, A = ±^^2^/3/15, and 
the typical frequency is / ~ 0.037. We will use these 
mean-field values later to compare with the simulation 
results. In the special case of symmetric reaction rates 
where ka = kb = kc = k, we get A = ±ifc/-\/3; for 
example, if /? = 1 and k ~ 0.5, then A — ±i^/3/6 and 
/ « 0.046. In addition, the system also has three absorb- 
ing states, with only a single species surviving ultimately: 
{p, 0, 0), (0, p, 0), and (0, 0, p). Within the mean-field ap- 
proximation, these fixed points are all linearly unstable. 
However, in any stochastic model realization on a finite 
lattice, temporal evolution would ultimately terminate 
in one of these absorbing states, as we shall explore for 
small systems below. 



III. MONTE CARLO SIMULATION RESULTS 

A. Model variants and quantities of interest 

We investigate stochastic RPS systems on one- and 
two-dimensional lattices with periodic boundary condi- 
tions. At each time step, one individual of any species 
is selected at random, then hops to a nearest-neighbor 
site, if the number of particles on the chosen target site 
is empty. Otherwise, one of the particles on the chosen 
neighboring site is selected randomly and undergoes a re- 
action with the center particle according to the scheme 
and rates specified by ^ if both particles are differ- 
ent. The outcome of the reaction then replaces the elim- 
inated particle. Note that predation reactions always in- 
volve neighboring particles; on-site reactions do not oc- 
cur. This has the advantage of permitting us to treat the 
model variants with and without site occupation num- 
ber restriction within the same setup, allowing for direct 
comparison. A similar approach was already adopted for 
the two-species lattice Lotka-Volterra model, where we 
confirmed earlier that nearest-neighbor predation inter- 
actions and strictly on-site reactions lead (without loss of 
gen erality) to essentially identical macroscopic features 

If the selected and focal particle are of the same species, 
the center particle just hops to its chosen neighboring 
site. For our model variants with site occupancy restric- 
tion, the hopping process only takes place if the total 
number of particles on the target site is less than the 
maximum occupancy number (local carrying capacity) 
Jirn- In this work, we set rim = 1; i.e., each lattice site 
can either be empty or occupied by a single particle of 
either species A, B,oi C (which gives four possible states 
for each site). Once on average each individual particle 
in the lattice has had the chance to react or move, one 
Monte Carlo step (MCS) is completed; thus the corre- 
sponding simulation time is increased by St ~ N~^. Also 
note that the hopping processes set the fundamental time 



Model 


Reaction rates 


Site restriction 


1 


homogeneous rate: k = 0.5 


no restriction 


2 


homogeneous rate: k = 0.5 


at most one particle 


3 


uniform rate distribution 


no restriction 


4 


uniform rate distribution 


at most one particle 



TABLE I: List of stochastic lattice RPS model variants. 



scale; basically the reaction rates are measured in units 
of the difFusivity D (unless D — 0). 

First, we shall study models with uniform symmetric 
reaction rates {ka = kb = kc = k = 0.5); next we sim- 
ulate systems with quenched spatial disorder by draw- 
ing the reaction probabilities k at each lattice site from 
a uniform distribution on the interval [0,1]. Therefore, 
this distribution has the same mean reaction rate 1/2 
as the homogeneous rate in the model with fixed reac- 
tion rates, allowing for direct comparison of the relevant 
numerical quantities. The four basic different model vari- 
ants we have investigated are summarized in Table HI 
In addition, we have studied systems with asymmetric 
reaction rates, both uniform and subject to quenched 
randomness with flat distribution. Besides the time- 
dependent population densities a(t), b{t), and c(i), aver- 
aged over typically 50 individual simulation runs, we also 
investigate their corresponding temporal Fourier trans- 
forms o(/) = / a{t) e^'^*^* dt, and the equal-time two- 
point occupation number correlation functions (cumu- 
lants) Cab{x, t) = {nA{i + x, t)nB{i, t)) — a{t) bit), where 
i denotes the site index, and similarly for the other 
species, as well as CAA{x,t), etc. In addition, for small 
systems with N lattice sites we have numerically com- 
puted the mean extinction time Teyi{N) defined as the av- 
erage time for the first of the three species to die out [1^ . 
For the one-dimensional four-state RPS model, we have 
also determined the time evolution of the typical single- 
species domain size (A(i)), see Sec. IIIIDl 

B. Two-dimensional stochastic RPS lattice models: 
symmetric rates 

We first report and discuss our Monte Carlo simula- 
tion results on a 256 x 256 square lattice with periodic 
boundary conditions. The data are typically averaged 
over 50 Monte Carlo runs with different initial config- 
urations, where the particles of each species are placed 
randomly on the lattice. Figure [Ta| depicts the tempo- 
ral evolution of the total population densities in a sys- 
tem without site occupation number restrictions and with 
equal reaction rates ka — kb — kc — 0.5 (labeled model 
1 in Table m, but unequal initial densities a(0) = 2/3, 
6(0) — c(0) = 1/6, along with two snapshots llbllcl of 
their spatial distribution at different times. Since the 
selection and reproduction processes are combined into 
a single step in our model, the total population den- 
sity p is strictly conserved, and as expected we there- 
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FIG. 1: (Color online.) (a) Temporal evolution for 
the population densities of species A (red/solid line), 
B (green/dashed), and C (blue/dash-dotted) with sym- 
metric reaction rates ka — kb — kc — 0.5 and without 
site occupation restriction (model 1), with unequal ini- 
tial densities a(0) = 2/3, &(0) = c(0) = 1/6, averaged 
over 50 Monte Carlo runs on a 256 x 256 square lattice, 
(b) Snapshot of the spatial particle distribution for a sin- 
gle simulation run at t = 50, and (c) a,tt = 500 MCS; (d): 
snapshot at t = 500 MCS, for a system where at most 
one particle of either species is allowed per site (model 
2). (Color coding shows the majority species on each 
site; red/gray: species A, yellow/light gray: B, blue/dark 
gray: C, black: empty site.) 



fore observe no spiral patterns that are characteristic of 
RPS models without conservation law [s^. In the initial 
time regime, we see distinct decaying population oscil- 
lations in Fig. Ila[ and inhomogeneous species clusters 
in the snapshot Fig. (Tb] As time progresses, the am- 
plitude of the oscillating fluctuations decreases quickly, 
and also the spatial distribution and species cluster size 
become more stable and homogeneous (Fig. [Tc)) . Our 
(fairly large) system eventually settles in a coexistence 
state with small density fluctuations (Fig. [Ta] inset). For 
comparison. Fig. Ildl shows a snapshot in a system with 
identical reaction rates and asymmetric initial densities, 
but with all site occupation numbers restricted to at most 
a single particle (model 2); one observes the same small 
cluster structure as in the absence of occupation restric- 
tions. 
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FIG. 2: (Color online). Signal Fourier transform |a(/)| 
of species A density data on a 256 x 256 square lattice 
with initial population densities a(0) = 6(0) = c(0) — 1/3 
for the four model variants described in Table HI averaged 
over 50 Monte Carlo simulation runs. 



In Fig. [5] we show the absolute values of the Fourier 
transformed population density signals |a(/)|, as ob- 
tained from averaging 50 Monte Carlo simulation runs for 
the four different model variants listed in Table HI in this 
case with equal initial densities a(0) = 6(0) = c(0) — 1/3. 
Recall that mean- field theory predicts a regular, un- 
damped oscillation frequency / ~ 0.046. From the sim- 
ulation data, we determine the characteristic peak fre- 
quency / « 0.028, which evidently governs oscillatory 
fluctuations; however, the finite width of the Fourier 
peak in Fig. [5] reflects that the population oscillations 
are damped and will cease after a finite characteristic 
relaxation time. 

Moreover, we see that even if spatial disorder and/or 
site occupancy restrictions are incorporated in the model, 
the Fourier-transformed density signals display practi- 
cally the same frequency distribution and significant peak 
locations. Indeed, we find that in our simulations for 
model versions 1 and 3 with total density 1, the typical 
occupation number at each site remains n < 2 through- 
out the runs, which explains why the exclusion con- 
straints in model variants 2 and 4 do not have a large 
effect. Thus, neither spatial disorder nor site occupancy 
restrictions change the temporal evolution pattern of the 
system markedly. This is in stark contrast with results for 
the two-species stochastic lattice Lotka-Volterra model, 
for which one finds (i) very pronounced spatio-temporal 
structures in the species coexistence regime [i^ ; (ii) large 
fluctuations that strongly renormalize the characteristic 
population oscillation frequency [l^, |4l| ; (iii) an extinc- 
tion threshold for the predator species induced by local 
density restrictions on the prey [40| ; and (iv) considerable 
enhancement of the asymptotic densities of both species 
caused by spatial variability of the predation rate . 

In order to study the effect of spatial disorder and 



5 




X X 



(a) (b) 

FIG. 3: (Color online.) (a) Static density autocorrela- 
tion function Caa{x) and (b) cross-correlation function 
Cab{x) (linear-Zogio plots) measured a.t t = 250 MCS for 
the four model variants described in Table HI with initial 
population densities a(0) = 6(0) = c(0) = 1/3. 



site occupation restriction on emerging correlations in 
our stochastic RPS models, we have determined the 
equal-time two-point correlation functions in the quasi- 
stationary (long-lived) coexistence state illustrated for 
models 1 and 2 in Figs. [Tc] and lldi respectively. These 
static correlation functions can quantitatively capture 
the emerging spatial structures in the lattice. Figures [3al 
and |3b] depict the autocorrelation function Caa (x) and 
the cross-correlation function Cab (x) as obtained for our 
four models (see Table H]), which all are seen to decay ex- 
ponentially with distance, i.e., CAAix) oc e"''^'/''*-'^ and 
Cab{x) oc e~l^l/'^^. From these log-normal plots, we 
have extracted the associated correlation length Iaa and 
typical species separation distance Iab] the results are 
listed in Table |lll It is worth noticing that in systems 
exhibiting spiralling patterns, as in the four-state RPS 
model without conservation law, the correlation functions 
Caa{x) and CABix) do not fall off exponentially but ex- 
hibit (damped) oscillations, see e.g. Refs. [sl,!!^. Site 
occupation restrictions clearly have the effect of reducing 
both correlation lengths. Also, as is the case for the two- 
species lattice Lotka-Volterra system [l^, rendering the 
reaction rate a quenched random variable for each site 
leads to more localized population and activity patches, 
characterized by markedly smaller correlation and typical 
separation lengths. 

The influence of varying (homogeneous and symmet- 
ric) reaction rates and modifying the total (conserved) 
population density is explored in Fig. |4l which shows the 
dependence of the characteristic Fourier peak frequency 
f on k and p. We find that / scales roughly linearly with 





Model 1 


Model 2 


Model 3 


Model 4 


Iaa 


3.27 ± 0.02 


2.92 ± 0.02 


2.64 ± 0.03 


2.59 ± 0.01 


Iab 


2.86 ± 0.08 


2.35 ± 0.09 


2.40 ± 0.06 


1.99 ± 0.09 



TABLE II: Correlation lengths Iaa for the autocorre- 
lation function and Iab for the cross-correlation func- 
tion (in units of the lattice spacing) obtained for the four 
model variants of Table U with symmetric reaction rates. 



0.05 I r 



0.04 - - 



0.03 - 




0.02 - ^ ' 

0.01 - a - ' ' ' 

0.00 I 1 ' 1 ' 1 ' 1 ' 1 ' 1 ' 1 

0.3 0.4 0.5 0.6 0.7 0.8 0.9 

k,P 

FIG. 4: ( Color online.) Variation of the characteristic 
peak frequency in the density Fourier signal |a(/)| with 
the total density p and homogeneous, symmetric reaction 
rate A:, for RPS simulations on a 256 x 256 square lattice 
with equal initial densities, run for 1000 MCS. 



both the total density p and the reaction rate A:, in ac- 
cord with the mean-field prediction f (x pk, see Sec. |lll 
We have also checked that switching off nearest-neighbor 
hopping (setting D — 0), thus allowing particle spreading 
only via the nonlinear reaction processes ([Ij, essentially 
leaves the stochastic RPS system's features intact. 

Finally, we have also studied the mean extinction time 
as function of lattice size N for small two-dimensional 
stochastic lattice RPS systems, here of the model 1 va- 
riety with homogeneous symmetric reaction rates ka = 
kb ^ kc — 0.5 and equal initial densities a(0) = &(0) = 
c(0) = 1/3. We recall that in any finite system display- 
ing an absorbing stationary state, stochastic fluctuations 
will eventually reach this absorbing configuration. In the 
stochastic RPS model, one therefore expects two species 
to eventually become extinct; however, reaching this ab- 
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(a) (b) 

FIG. 5: (a) Mean extinction time as function of lattice 
size N (linear-iogio plot), obtained from averages over 
50 Monte Carlo runs, for small two-dimensional lattice 
RPS systems in the absence of site restrictions and with 
symmetric reaction rates ka — kb — kc — 0.5 (model 
1), and equal initial population densities a(0) — 6(0) = 
c(0) = 1/3. The data are for lattices with = 5 x 5, 
7 X 7, 10 X 10, 12 X 12, 15 x 15, 17 x 17, and 20 x 20 sites, 
(b) Histogram of measured extinction times for A^ = 100 

sites. 
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sorbing state may take an enormous amount of time, 
and will thus become practially unobservable on large 
lattices. In fact, in two and higher dimensions one ex- 
pects the mean extinction time Tox (here measured for 
the first species to die out) to scale exponentially with 
system size N, since random fluctuations effectively have 
to overcome a finite barrier in order to follow an 'opti- 
mal' path towards extinction. As depicted in Fig. I5a[ 
we indeed observe \nTox{N) ^ N, consistent with the 
prediction on the coexistence state stability reported in 
Refs. [H, [3^. The associated distributions of extinction 
times are described by neither Poisson nor Gaussian dis- 
tributions (e.g., the means are considerably larger than 
the most likely values), but display long 'fat' tails at large 
extinction times, see Fig. I5bl We expect similar features 
in model variant 2, in accord with the remarkably long- 
live species coexistence observed in Rcf. [23 |. 



C. Two-dimensional stochastic RPS system: 
asymmetric rates 

Next we turn to a stochastic RPS system with asym- 
metric reaction rates and consider the various model vari- 
ants specified in Table HIT] together with the reactions ([ij. 
Figure I5al shows the time evolution for the three species' 
densities in a system with constant rates ka = 0.2, kb ~ 

0. 5, and fee = 0.8. From our simulations for model version 

1, we infer the asymptotic population densities (with sta- 
tistical errors) (0.40±0.01, 0.45±0.01, 0.15±0.01), which 
follow the trends of the mean-field results (a* ,b*,c*) — 
(0.33,0.53,0.13). As becomes apparent in the snapshots 
I6bl and [6c] for model variant 1 without site restrictions, 
and I6dl for a system with at most a single particle per 
site (model 2), particles of the same species form dis- 
tinctive spatial clusters. The effect of the reaction rate 
asymmetry on the equal-time auto- and cross-correlation 
functions is shown in Figs. [7al and I7bl respectively, with 
the ensuing correlation lengths and typical separation 
distances listed in Table IIVI Note that the autocorre- 
lation length Icc for species C is smaller than Iaa, and 
Ibb, which is largest. This is consistent with the long- 
time densities in the (quasi-stationary) coexistence state, 
given our observation that the overall particle density is 



Model 


Reaction rates 


Site restriction 


1 


ka = 0.2,^6 = 0.5, kc = 0.8 


no restriction 


2 


ka = 0.2,^6 = 0.5, kc = 0.8 


at most one 


3 


ka e [0,0.4], fcfc = 0.5, kc = 0.8 


no restriction 


4 


ka e [0,0.4], fci, = 0.5, kc = 0.8 


at most one 



TABLE III: List of stochastic lattice RPS model variants 
with asymmetric rates. While kf, = 0.5 and kc = 0.8 are 
held fixed in all four variants, we set ka — 0.2 in models 
1 and 2, whereas we took ka to be a random variable 
uniformly distributed in [0, 0.4] in models 3 and 4. 
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FIG. 6: (Color online.) (a) Temporal evolution for 
the population densities of species A (red/solid line), B 
(green/dashed), and C (blue/dash-dotted) with asym- 
metric reaction rates ka — 0.2, kb = 0.5, kc = 0.8 and 
without site occupation restriction (model 1), with equal 
initial densities a(0) — 6(0) = c(0) = 1/3, averaged over 
50 runs on a 256 x 256 square lattice, (b) Snapshot of the 
spatial particle distribution in a single simulation run at 
t = 50, and (c) at t = 500 MCS; (d): snapshot at t = 500 
MGS, for a system where at most one particle of either 
species is allowed per site (model 2). (Majority species 
coloring: red/gray: A, yellow/light gray: B, blue/dark 
gray: C, black: empty.) 



roughly uniform. 

As a last model variation, we allow the reaction rate ka 




FIG. 7: ( Color online.) (a) Equal-time autocorrelation 
functions Caa{x), Cbb{x), Ccc{x) at t = 1000 MGS 
for the model described in Fig. 6. (b) Equal-time cross- 
correlation functions CAsix), Cbc{x), Cac{x) (linear- 
^0510 plots). 
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Iaa Ibb Icc 

5.24 ± 0.03 5.67 ± 0.08 3.46±0.05 

Iab Ibc Iac 

6.68 ± 0.20 3.68 ± 0.07 3.33 ± 0.05 



TABLE IV: Correlation lengths (top) inferred from 
the autocorrelation functions and typical separation dis- 
tances (bottom) obtained from the cross-correlation func- 
tions (in units of the lattice spacing) measured for the 
RPS model with asymmetric but homogeneous reaction 
rates ka = 0.2, fc(, = 0.5, and kc = 0.8. 



0.8 




200 400 600 800 1000 

t 




(a) (b) 

FIG. 10: (Color online.) (a) Equal-time autocorrela- 
tion function Caa (x) and (b) cross-correlation functions 
Cab{x) (iineai-logio plots) at t = 1000 MCS for the four 
model variants described in Table Hill 

but hold kb = 0.5 and kc = 0.8 fixed. Fig. |S] compares 
the time evolution for these disordered systems with and 
without site restrictions with the corresponding homoge- 
neous models. Once again, we see that spatial variability 
in the reaction rate even in this asymmetric setting has 
very little effect. As can be seen from the Fourier signal 
peak in Fig. ^ the characteristic frequency comes out to 
be / « 0.021 for all four asymmetric model variants in- 
vestigated here, and Figs. llOal and llObl demonstrate that 
the disorder hardly modifies the spatial decay of the auto- 
and cross-correlation functions either. 



FIG. 8: (Color online.) Time evolution for the popu- 
lation density a{t) of species A for four model variants 
with asymmetric reaction rates, namely with k^ = 0.5, 
kc ~ 0.8 and either uniformly ka = 0.2, or drawn from 
a flat distribution [0,0.4], with and without site restric- 
tions (see the listing in Table HIT)) . The initial densities 
are a(0) = 6(0) = c(0) = 1/3, and the data stem from 
averages over 50 runs on a 256 x 256 square lattice. 




FIG. 9: (Color online.) Signal Fourier transform |a(/)| 
for the four model variants described in Table IIIIl The 
characteristic frequency comes out to be / w 0.021. 



to be a quenched spatial random variable drawn from the 
flat distribution [0,0.4], such that its average is still 0.2, 



D. One-dimensional Monte Carlo simulations 

We have run simulations for all four model variants 
listed in Table U i.e., with/ without site occupancy re- 
striction; with/without quenched spatial randomness in 
the reaction rates, in one dimension. We find that only a 
single species ultimately survives and eventually occupies 
the whole lattice no matter whether spatial disorder or 
site restrictions are included in the model: as expected, 
the one-dimensional system will always evolve towards 
one of the three extinction states where two of the three 
species will die out. While this phenomenon also occurs 
in two dimensions, in d = 1 the species coexist over a 
time that on average scales polynomially with the sys- 
tem size (see below), i.e., extinction happens on a much 
shorter time scale than in two dimensions (see Fig. [Sa]) . 
Again, for equal (mean) reaction rates and initial den- 
sities, each species has equal survival probability. For 
comparison, the space-time plots of one-dimensional lat- 
tice simulations without and with site occupancy restric- 
tion are depicted in Figs. [TT] and 12, respectively. It is 
seen that individuals of identical species cluster together, 
and any reactions are confined to the boundary separat- 
ing the single-species domains. When the occupancy of 
any site is restricted to a single particle of either species, 
these domains form quickly and are very robust, even if 
not all sites are filled, see Figs. I12al and [T2bl 

The population density signal Fourier transform [a(/)[, 
shown for species A in Fig. I13[ confirms the absence of 
any population oscillations through the absence of any 
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FIG. 11: (Color online.) Time evolution (up to 1000 
Monte Carlo steps; from top to bottom) for a one- 
dimensional RPS model run with equal, homogeneous 
reaction rates ka = kh = kc = 0.5, equal initial densities 
a(0) = 5(0) = c(0) = 1/3, and in the absence of site occu- 
pancy restriction. (Only 10000 of the total 50000 lattice 
sites in this run are shown; majority species coloring in 
red/gray: A, yellow/light gray: B, blue/dark gray: C, 
black: empty.) 




0.010 ^ 0.015 



FIG. 13: (Color online.) Signal Fourier transform |a(/)| 
for the four RPS model variants listed in Table HI in one 
dimension. The data is averaged over 50 Monte Carlo 
simulations on lattices with 50000 sites. 




FIG. 12: (Color online.) Time evolution (up to 
1000 Monte Carlo steps; from top to bottom) for one- 
dimensional RPS model runs with equal, homogeneous 
reaction rates ka = k^ = kc ~ 0.5, equal initial densities 
(a) model 2: a(0) = 5(0) = c(0) = 1/3, and (b) model 
2': a(0) = 5(0) = c(0) = 0.2 (model 2' and 4' refer to the 
corresponding model variants listed in Table H] with total 
particle density less than 1), where at most one particle of 
either species is allowed per site. (Only 10000 of the total 
50000 lattice sites in these runs are shown; red/gray: A, 
yellow/light gray: B, blue/dark gray: C, black: empty.) 



Ref. [2J]. The domain stability is further illustrated by 
the very slow temporal decay of the on-site auto- and 
cross-correlation functions fsee Figs. I15al and I15bp . No- 
tice that quenched spatial disorder in the reaction rates 
does not affect the time evolution of the autocorrelation 
functions, in contrast with site occupancy restrictions; 
here the results depend on the presence or absence of 
empty sites, see Fig. I15al However, the cross-correlation 
functions in Fig. ll5bl look essentially indistinguishable for 
all these model variations. 



A 
V 



- with site occupancy restriction 
withiout site occupancy restriction 



peak at nonzero frequency /, and the width of the peak 
at / = reflects the decay time to the stationary ex- 
tinction state. As in two dimensions, we observe very 
little effect of either site occupation number restrictions 
or spatial variability of the reaction rates on the Fourier 
signal, compare Figs. [2] and [HI We have also measured 
the mean single-species domain size (A(t)) and investi- 
gated its growth with time t, shown in Fig. [TH As was 
predicted in Refs. [25, 26], for the implementation with 
site occupancy restriction (model 2) to at most a sin- 
gle particle per site, we observe (A(t)) ~ t^^^] we find 
the same asymptotic growth law when arbitrarily many 
particles are allowed on each lattice site. An algebraic 
decay of the number of domains was also reported in 



1 1 ' — ' — 1 

10 100 1000 

t 

FIG. 14: ( Color online.) The time evolution {logiQ—logiQ 
plot) of the mean single-species domain size (A(i)) mea- 
sured in one-dimensional Monte Carlo simulation runs 
with 10000 lattice sites for RPS models with symmetric 
reaction rates ka — kh = kc = 0.5 and equal initial den- 
sities a(0) — 5(0) = c(0) = 1/3. The upper (black) curve 
shows the data for model 2 with site occupation restric- 
tion, see Fig. I12al whereas the lower (red/dashed) graph 
pertains to model variant 1 without site occupancy re- 
strictions, see Fig. [TTl For comparison, the blue/dotted 
straight line represents the predicted t^^^ power law. 



(a) (b) 

FIG. 15: (Color online.) Time evolution for the on-site 
(a) autocorrelation CAA{0,t) and (b) cross-correlation 
function CAB{Q,t) in one-dimensional RPS model vari- 
ants with 500 sites, averaged over 1000 simulation runs. 
Shown are the results for model variants 1-4 with ini- 
tial densities a(0) ~ 6(0) = c(0) = 1/3; model 2' refers 
to a system with site occupancy restriction 1 and lower 
particle density a(0) = b{0) = c(0) = 1/4. 



Figures I16al and I16bl respectively depict the equal- 
time auto- and cross-correlation functions for the vari- 
ous model variants listed in Table |T] obtained for a one- 
dimensional lattice with 50000 sites. We observe expo- 
nential decay with similar large correlation lengths for all 
model variants upto about 50 lattice sites, followed by a 
cutoff (which extends to larger x as time increases). 

Finally, we investigate the mean extinction time as 
function of system size N in one dimension. As becomes 
apparent in Figs. I17a[ in all one-dimensional model vari- 
ants we have considered, within our (large) error bars 
the mean extinction time appears to follow a power law 
Tex - iV^, as proposed in Refs. [H [H S, . However, 
a best power-law fit yields variable effective exponents 
Tex ~ N'^ with 7 ranging from ~ 1.5 to ~ 1.8 if we fit 
the data up to = 50 or N = 200, respectively, rather 
than 7 = 2 or 7 = 1 [lH, HI, Egl . Biasing the data 
towards smaller systems for which the statistical errors 
are likely better controlled, our results may even be con- 
sistent with the mean- field value 7 = 1. Note, however, 
that the extinction time distribution acquires even fatter 




™ "° "° ° x 

(a) (b) 



FIG. 16: (Color online.) (a) Static autocorrelation func- 
tions Caa{x) and (b) static cross-correlation functions 
Cab{x) (linear-/o(7io plots) measured at t = 250 MCS 
for the four RPS model variants described in Table U on 
a one-dimensional lattice with 50000 sites, averaged over 
50 simulation runs. 



(a) (b) 

FIG. 17: (Color online.) (a) Mean extinction time as 
function of lattice size N {logio — logio plot), obtained 
from averages over 1000 Monte Carlo runs, in one dimen- 
sion, (b) Histogram of the measured extinction times for 
model variant 2 with N = 30; compare with Figs. l5al and 

m 



tails at large times than in two dimensions, see Fig. I17bl 
and rare long survival events dominate the averages and 
induce large statistical fiuctuations. The mean extinction 
time alone therefore poorly characterizes the extinction 
kinetics. When the reaction rates are chosen asymmet- 
ric, we have checked that only the 'weakest' species with 
the smallest predation rate survives^whereas the other 
two species are driven to extinction [23, 1231 ■ 



IV. CONCLUSION 

In this paper, we have studied the effects of finite car- 
rying capacity and spatial variability in the reaction rates 
on the dynamics of a class of spatial rock-paper-scissors 
(RPS) models. We have investigated the properties of 
several variants of the stochastic four-state zero-sum RPS 
game (with conserved total particle number) on two- and 
one-dimensional lattices with periodic boundary condi- 
tions. In two dimensions, owing to the strict (local) con- 
servation of the total particle number, one does not ob- 
serve the formation of spiral patterns; the three species 
simply form small clusters. In fact, spatial correlations 
are weak in the (quasi-stationary) coexistence state, and 
the system is remarkably well described by the mean-field 
rate equation approximation. Typical extinction times 
scale exponentially with system size (36j , resulting in co- 
existence of all three species already on moderately large 
lattices. We find the characteristic initial oscillation fre- 
quency to be proportional to the reaction rate and total 
particle density, as predicted by mean-field theory. 

We observe that neither site occupation number re- 
strictions nor quenched spatial disorder in the reaction 
rates markedly modify the populations' temporal evolu- 
tion, species density Fourier signals, or equal-time spa- 
tial correlation functions. This observation holds for 
models with symmetric as well as asymmetric reaction 
rates, and even if spatial variability is introduced only 
for the competition of one species pair. On the basis 
of the mean-field results, this very weak disorder effect 
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is a consequence of the essentially linear dependence of 
the long-time densities on the reaction rates fc; averag- 
ing over a symmetric distribution just yields the average. 
In the two-species Lotka-Volterra model, instead both 
the asymptotic predator and prey densities are inversely 
proportional to the predation rate, and averaging over a 
distribution of the latter strongly biases towards small 
rate values and large densities [42,]. Both in the two- 
and cyclic three-species systems, spatially variable rates 
induce stronger localization of the species clusters. 

In one dimension, two species are driven towards ex- 
tinction with the mean extinction time Tox ~ N'^ , 7 « 
1 ... 1.8 (with large error bars), and only a single species 
survives. The distribution of extinction times displays 
fat long-time tails. We confirm that the single-species 
domains grow with the predicted power law {X{t)) ~ t^^^ 
in the models with site restrictions [2^, [2^ , and we have 
obtained similar results for the model variants with an in- 
finite carrying capacity. For asymmetric reaction rates, 
we have also checked that the 'weakest' species is the 
surviving one [H, [2^ . 

Our results demonstrate that the physical properties of 
cyclic RPS models are quite robust, even quantitatively, 
with respect to modifications of their 'microscopic' model 
definitions and characterization. This is in stark contrast 
with the related two-species Lotka-Volterra predator- 
prey interaction model. We believe the origin of this 



remarkable robustness lies in the comparatively weaker 
prominence of stochastic fluctuations and spatial correla- 
tions in the cyclic three-species system. The robustness 
of the RPS models considered here implies that environ- 
mental noise can be safely ignored and that their prop- 
erties are essentially independent of the carrying capac- 
ity. It is worth emphasizing that this result is nontriv- 
ial; in terms of modeling such systems, it notably implies 
that one has the freedom to consider strict site restriction 
("■m = 1) E^nd hence simplify the numerical calculations, 
or to set n„ — 00 (no site restrictions) and thus facili- 
tate the mathematical treatment. In fact, our study es- 
tablishes that both 'microscopic' model realizations are 
essentially equivalent. As it turns out, this conclusion 
also pertains to spatial May-Leonard models [i^ 113] , as 
we shall report in detail elsewhere [45j . 
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